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SOCIALLY OPTIMAL CHARGING STRATEGIES FOR 
ELECTRIC VEHICLES 
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. Abstract. Electric vehicles represent a promising technology for 

Q^ reducing emissions and dependence on fossil fuels and have started 

•^T entering different automotive markets. In order to bolster their 

adoption by consumers and hence enhance their penetration rate, 
a charging station infrastructure needs to be deployed. This pa- 
i— i per studies decentralized policies that assign electric vehicles to a 

Q.J network of charging stations with the goal to achieve little to no 

^^ queueing. This objective is especially important for electric vehi- 

cles, whose charging times are fairly long. The social optimality 
of the proposed policies is established in the many-server regime, 
where each station is equipped with multiple charging slots. Fur- 
ther, convergence issues of the algorithm that achieves the optimal 
policy are examined. Finally, the results provide insight on how 
,__! to address questions related to the optimal location deployment of 

^- the infrastructure. 

On 

CO 1. Introduction and model 

(N 

.• There has been an increasing penetration of Plug-in Hybrid and pure 

Q Electric Vehicles (PHEV/EV) over the least few years pQ . This is due 

to developments in battery technology that have dramatically increased 
their range [2], advances in charging technologies that have reduced 
their charging times pJ], incentives that have lowered their acquisition 

/\ and operation cost and an overall desire to lower emissions [I]. 

On the other hand, they represent a potential source of disruption 
to normal grid operations if not integrated carefully because they need 
to connect to the distribution network to charge [HJ|6]. Currently, EVs 
are equipped primarily with lithium-ion batteries, ranging in energy 
capacity from 5 kWh for short-range PHEVs to 50 kWh for high per- 
formance EVs. Further, today's and future EVs are designed with a 
wide range of specifications to satisfy different customer preferences. 
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Their impact on power grid operations will heavily depend on their 
market penetration [3]. Estimates vary widely, ranging from 3 to 18 
million vehicles by 2025 and from 5 all the way to 40 million vehicles 
(approximately 20% of the total US market) by 2030 [4|. Their dis- 
ruptive impact is mainly due to the energy load they represent. On 
average, under normal charging conditions (1.4 kW) an EV represents 
a 1.3-fold of a full household load, whereas fast charging technologies 
(7.2 kW) correspond to an almost 3-fold increase (HE]. The injection 
of such large loads coupled with the possible uneven geographic dis- 
tribution of EVs would definitely strain the entire grid as argued in 

Given the current predominance of PHEV vehicles, the literature has 
largely focused on scheduling at home overnight charging (see [HJ [10] 
and references therein). The proposed approaches treat the induced 
load as an aggregate and discuss different mechanisms on how to shift 
it during night hours to take advantage of the underutilized electricity 
generation assets. However, with increasing penetration rates, efficient 
operation of an expanding charging station infrastructure becomes a 
key issue. Charging EVs is a rather slow process, as even fast charging 
takes at least half an hour [4], thus requiring careful scheduling poli- 
cies to provide the necessary quality of service to customers. Faster 
charging technologies (e.g. DC charging) could mitigate some of these 
effects, but as mentioned above, electric utilities have concerns about 
possible negative impacts of such technologies on the power grid, if 
deployed at large scale. The work to date on charging stations has 
mostly focused on modeling and optimizing the architecture of a single 
charging station (TTJ [12], [13] . 

In this paper, we consider decentralized dispatching policies that as- 
sign electric vehicles to a network of charging stations. Our focus is 
on guaranteeing quality-of-service to charging customers. Due to the 
lengthy charging times, our primary focus is to ensure little to no cus- 
tomer queueing. Our results apply to the "many-server" regime, when 
each charging station is equipped with multiple charging slots (e.g. 
a commercial parking lot, a strip of street parking in a densely urban 
area), possibly in a mix of charging technologies (e.g. Level-2 charging, 
coupled with very fast DC charging infrastructure). 

We provide next a description of the modeling framework. Consider 
EVs and a network of charging stations within a specific geographical 
region, for example an urban area or a section of the highway system. A 
vehicle that needs to recharge broadcasts a signal indicating its location 
and battery type and status, and receives responses from the charging 
stations in its neighborhood. We assume that the EV has preferences 
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among the different charging stations, which we model through costs; 
these may be related to the distance of the current location of the EV 
to the charging station, or simply whether it is capable of reaching the 
station using its remaining battery power. Based on the stations' re- 
sponses, the EV chooses a charging station and immediately proceeds 
there. We are interested in designing socially optimal charging strate- 
gies; namely, to engineer the response signals from stations to vehicles 
and the resulting routing decisions, so that the average cost incurred 
by the vehicles is minimized. 

We will make the simplifying assumption that the vehicles can be 
partitioned into a finite set of types % = 1, . . . ,1 which encode their 
location, preferences, and battery technology. Thus, the signal that 
the EV broadcasts to the nearby charging stations indicates the type 
of the vehicle. We assume also that the various charging slots available 
in the network can be partitioned into a finite set of types j — 1, . . . , J, 
encoding their location and technology. We are interested in the regime 
where there are many charging slots of each type: say Nj chargers of 
type j, with Nj large (e.g. at least 10). Thus, as mentioned above, 
the architecture of the charging station in the network is such that it 
comprises of many identical chargers. 

Let Xi be the rate at which EVs of type % make charging requests. 
Let /i^- 1 be the expected time it takes for a charger of type j to satisfy a 
request of type i; this may be infinite (corresponding to Hij = 0), e.g. if 
an EV cannot reach the charger or their technologies are incompatible. 
Finally, let q(j) be the cost that an EV of type % incurs on being 
assigned to type j; this may be a measure of the distance between the 
vehicle and the chargers of type j, or a constant whenever the vehicle 
can reach the charger. We set q(j') = oo whenever an EV of type 
% can not reach a charger of type j. Without any guidance from the 
system/network, we would expect a vehicle of type i to always choose a 
charger of type j that minimizes q(j); however, the following example 
shows this behavior to be suboptimal in some cases. 

Consider a system with two vehicle classes A and B, and two charger 
classes 1 and 2. Suppose that A cannot reach class 2, and B can use 
both charger classes, but prefers type 1. That is, ca(2) = oo and 
cb(1) < cb(2). In this case, if both EV types preferentially go to 
chargers of type 1, vehicles of type A may have to queue. If, on the 
other hand, vehicles of type B can be directed to chargers of type 2 
when few chargers of type 1 remain, queueing can be avoided. 

Another example is obtained by considering two EV and charger 
classes as before. Suppose that A prefers type 1 and B prefers type 
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2; that is, c^(l) < ca(2) and cb(2) < cb(1). Suppose further that 
/Mi < HA2 while /isi = hb2, meaning that it is faster to charge vehicles 
of type A at the more distant charger. Clearly, in heavy traffic we 
must encourage vehicles of type A to charge at the more distant, fast 
charging station 2, while vehicles of type B will correspondingly need 
to be directed to the more distant (and not even faster!) charging 
station 1. 

From a social perspective, the main objective is to to design charging 
strategies for a network of stations, where such routing decisions will 
happen automatically. 

This problem has similarities to the inventory and facility location 
problem: we are interested in "distributing" a finite supply of chargers 
among vehicles, subject to location constraints, in such a way as to 
keep vehicles from being "undersupplied" . Extensive literature exists 
on such problems; see, for example, [H] and references therein. How- 
ever, there are two key differences in our set-up. The first is that EV 
demand is mobile. Thus, demand from the vehicles of type (location) 
i can be split between several different charging stations (facilities). In 
fact, our work shows that this splitting should be encouraged. The 
implication from an algorithmic perspective is that instead of solving 
an integer program that typically the case in the standard facility lo- 
cation problem, we need to solve its convex relaxation. The second 
key difference between our setting and the inventory / facility location 
problem is that we wish to avoid centralized decisions. Instead, for our 
set-up to be scalable to the size of a large urban or even a metropolitan 
area with hundreds of charging stations and thousands of vehicles, we 
need the decision-making to be distributed. The latter goal justifies our 
formulation for letting the EVs choose their preferred charging station 
based on the information communicated to them by nearby stations, 
rather than involving some centralized planning scheme. 

To achieve this goal and design the required efficient, distributed al- 
gorithm for routing EVs to charging stations so as to avoid excessive 
delays due to queueing, we employ ideas from queueing theory and 
communication systems. We introduce the GPD algorithm (Greedy 
Primal-Dual) that has been successfully used in the call center litera- 
ture and establish analytic guarantees for its performance in large scale 
network of charging stations. Its main feature is its distributed an on- 
line nature, and also its automatic adaptability to changing arrival 
patterns. We also present two variants (Load Balancing and Freest 
Charger Shortest Queue), which are shown to exhibit superior perfor- 
mance to the GPD algorithm in selected settings. 
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1.1. Modeling Assumptions. We discuss next the main modeling 
assumptions. First, we have stochastic assumptions on key processes 
of the system under consideration, that lead to the analytic guaran- 
tees on the system behavior. Second, we specify which of the system 
parameters are known to which participants; the latter assumptions 
ensure that the system behaves in a distributed manner. 
We begin with stochastic process assumptions. 

Assumption 1. The arrival process of charging requests of EVs of 
type i satisfies a functional law of large numbers approximation. Let 
Ai(f) = ^{requests of type i up to time t}; we require 

-AAri) ==>■ Ajt, as r — >■ oo, uniformly on compact sets, 
r 

Throughout, the notation ==>■ incicates uniform convergence on com- 
pact sets. Furthermore, we require that arrivals have bounded second 
moment: K[(Ai(t + 1) — Ai(t)) ] < oo uniformly in t. 

The process of service completions of EVs of type i by chargers of 
type j satisfies a functional law of large numbers approximation. Specif- 
ically, let Sij(t) be the number of EVs of type i that have completed 
service with one charger of type j when that charger has spent a total 
amount of time t charging vehicles of type i; we require 

-Sij{rt) =>■ fiyt, as r — >■ oo ; u.o.c. 

We also require that service times have bounded second moment: K[(Sij(t+ 
1) — Sij(t)) } < oo uniformly in t. 

These are quite general assumptions, that are satisfied for example 
when interarrival and service times are independent and identically 
distributed (iid) and possessing finite second moments. 

Assumption [T] will be taken to hold throughout the paper. It is suf- 
ficient to show that the proposed algorithm has optimal throughput, 
meaning that it will successfully recharge all EVs whenever it is pos- 
sible to do so (perhaps with queueing delays). However, in order to 
provide more precise guarantees, for example on the probability of an 
arriving vehicle finding a free charger, we will need more control over 
the deviations of the arrival and service processes from their FLLN 
approximation. Assumption [2] asserts that the arrival and service pro- 
cesses obey a functional central limit theorem. 
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Assumption 2. The arrival processes of requests of type i obey the 
functional central limit theorem: 

:(Ai(rt) — XiTCJ =>■ Wi(t), as r — >• oo, u.o.c. 



Here, Wi are a set of independent Brownian motions (one per EV type) 
with some finite variance. 

In addition, the service processes of EVs of type i by chargers of type 
j obey the functional central limit theorem: 

: (Sij (rt) — fiijrt) =>■ Wij(t), as r — >■ oo, u.o.c. 



Here, Wij are a set of independent Brownian motions (one per pair of 
vehicle type and charger type) with some finite variance. 

These assumptions are standard in the queueing literature. They 
hold when the interarrival and service times are iid with finite sec- 
ond moment, and typically allow one to model the queueing process 
as a reflected Brownian motion. We will comment further on this in 
Section 13.31 

We next specify what information is available to which participants. 
The system has parameters A» (arrival rate of requests of type i), Hij 
(service rate of requests of type i by chargers of type j), Nj (number 
of chargers of type j), and q(j) (the cost associated with assigning a 
vehicle of type i to a charger of type j). In addition, the scheduling 
algorithm will use a parameter f3. 

Assumption 3. The number of possible types is finite. 

The parameters A« are unknown to anyone in the system; the algorithm 

will implicitly estimate them. 

The parameters fiij are assumed to be known to a request for which they 

are relevant (e.g., included as part of the exchange between the EV and 

the charging stations); if /x^- is not communicated, it is assumed to be 

0. 

The costs Ci(j) are assumed to be known to the EV for which it is 

relevant. Specifically, we require the EV to be able to compare quantities 

of the form c%{j) + KijUjj , where K^ will be a quantity communicated 

by the charging station. 

The parameter (3 (a small real number) is assumed to be the same at 

all charging stations. 

Note that all communications and knowledge are local: EVs only 
acquire information about nearby stations, and stations only acquire 
information about the nearby requests. 
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The remainder of the paper is organized as follows: In Section [2], we 
formulate the main algorithm, GPD, and show its throughput optimal- 
ly. In Section [3j we present a detailed analysis of the behavior of large 
systems, provide Brownian motion approximations for queue sizes and 
the associated probability of finding a free charger upon arrival. We 
also state a key result (Theorem [TJ) whose proof is given in the Appen- 
dix. In Section El we introduce the Load Balancing (LB) algorithm, 
which is designed to spread the load more evenly between chargers to 
avoid excess queueing. Section [5] discusses the case when the costs 
Ci(j) G {0, oo}, i.e. EVs are indifferent between chargers provided they 
are within reach. For the case of no user costs, we present another 
algorithm, Freest Charger Shortest Queue (FCSQ), which reacts faster 
than GPD or LB to changes in the arrival pattern. Section [6] presents 
selected simulation results of the behavior of the three algorithms on a 
simple system. Finally, Section [7] discussed insights into planning the 
charging network, as well as directions for future research. 



2. Scheduling policies and throughput optimality 

In this section, we formulate the policy that will achieve socially op- 
timal average costs while keeping queueing low, if possible. We begin 
by formulating the corresponding linear program. Let Ay be the av- 
erage rate at which vehicles of type i are routed to stations of type j. 
The basic linear program we consider is as follows. 

minimize /"AyC^") (l a ) 

X i = y %2\ ij ,Vi (lb) 



s.t. 



Nj > E — ' *J (1c) 

i Vij 

over Ay > 0, Vi, j. (Id) 

The objective function here is the rate at which costs are incurred. The 
constraints are the basic feasibility ones: all arriving requests need to 
be assigned, but on average no more than Nj chargers of type j may be 
used. Note that the optimal solution may be infinite, corresponding to 
insufficient capacity in the system. Call the value of this linear program 
S(X). 

We define the feasible region to be the set of routing rates for which 
the solution to this linear program is finite and the capacity constraints 
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(lc) are strictly satisfied. That is, 



A = {(k)i=i,...,i : S(X) < oo, and Nj > J^ ~ ^ 

j Hi 

for some optimal solution X*j of ([I]). 

We now introduce the Greedy Primal-Dual (GPD) Algorithm of 
which will determine a scheduling rule that implicitly solves the linear 
program (fij). The parameter f3 is a small real number that deter- 
mines the trade-off between convergence speed and precision of the 
solution. While the algorithm is a special case of the one in [15] . tech- 
nical differences between our setting and that in [15] imply that analytic 



guarantees presented in Section |3.3| are not automatic and need to be 
established rigorously. 



GPD Algorithm: 

1. Each charger type maintains a virtual queue variable Qj(t) (this 
need not be an integer). The latter are appropriately initialized, 
e.g. Qj(0) = 0. 

2. When an EV of type % requests service at time t, we locate 



G argmin Cj(j) + 



PQi(t) 



3 Hij 

That is, all stations in the neighborhood of i communicate pa- 
rameters /3Qj(t) (or (3Qj(t)/fiij), and the user picks the station 
j* that minimizes Ci(j) + f3Qj(t)/fiij. The vehicle announces 
its decision to j*, and the corresponding virtual queue is incre- 
mented: 



3. Decrease all virtual queues at rate Nj per time unit, provided 
they are positive. (Once a virtual queue hits 0, it stays at level 
until some EV is routed to it.) 

Note that the algorithm runs in continuous time and no synchronization 
is necessary between different stations. 

One important instance of the algorithm is when Ci(j) G {0, oo} for 
all i,j; that is, the EVs are indifferent among the charging stations as 
long as they can use them at all. We will discuss this case in detail in 
Section |5j For now we point out that in this case, the parameter (3 is 
unnecessary. 
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The results of [15] imply that whenever A ^0, the GPD algorithm 
will stabilize all virtual queues. This intuitively means that the algo- 
rithm will be routing vehicles to chargers at rates Ay G A. The more 
precise statement is as follows. Let Ay(£) denote the number of vehi- 
cles of type i that have been routed to chargers of type j up to time 
t. The limit lim^oo ^Aij(t), denoting the average rate, may not exist; 
however, any sequence of times t — > oo has a subsequence along which 
such the limit does exist. Pick one such set of subsequential limits, 
and call it (Ay). Then, if A ^ 0, the rates Ay satisfy the routing 



constraint (lb) and the capacity constraint (lc). In general, we may 
have J2ij ^ij c i(j) > S(X), i.e. the rates may not be socially optimal; 
note, however, that the algorithm remains throughput-optimal, that 
is, it stabilizes the system whenever the arrival rates belong to A. In 
Section [3] we will see that as (3 — > 0, the algorithm converges to the 
socially optimal routing rates. 

Note that in the case of Poisson arrival and service process, and with 
service rates /^y being rational numbers, the virtual queueing system 
is a countable state-space Markov process. In that case, there will 
be well-defined steady-state rates Ay to which the above observations 
apply. 

Further, note the natural form (from the user's point of view) of the 
routing decision: the EV driver is asked to add to her intrinsic costs 
Ci(j) a certain charge per unit time (3Qj(t). This charge will be greater 
for the stations that are in high demand, and lower for the stations that 
are less congested. We point out that, because the GPD algorithm uses 
virtual queues rather than actual queues to make routing decisions, it 
can be run in the background of some other scheduling mechanism, 
provided a model for the costs q(-) is available. 

3. Limiting regimes 

We now describe the asymptotic behavior of the system in certain 
limiting regimes. Our interest is in large systems; consequently, we will 
be interested in the case Aj — > oo. To accommodate the increasing 
arrival rate, we will consider Nj — > oo (many chargers at each sta- 
tion), holding the number of charger types and the speed of charging 
fixed. We will also consider the effect of taking /? — > 0, where is the 
parameter used in defining the GPD algorithm. 

In order to state the results, we make an assumption on the solution 
structure of the linear program (fTl) . Some of the conclusions apply even 
when these do not hold; however, the exposition would become more 
cumbersome. 
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Assumption 4. The optimal solution (A*) to the linear program is 
unique. The optimal set of dual variables (q*) is unique. The pairs (ij) 
for which A*,- > are called basic activities. 

This is a part of the complete resource pooling assumption, which 
is commonly made in many-server queueing literature, e.g. [16J. As 
pointed out in [16], Assumption [I] holds for generic values of parameters 
X i:j , fi i:j , and Ci(j). 

We describe briefly the techniques used to obtain the stability and 



"fluid-scaled" convergence results in Sections 3.1-3.2, relying only on 
Assumption [l] for the underlying stochastic processes. The technique is 
standard in the queueing literature, and involves the use of fluid limits 
rather than constructions of explicit Lyapunov functions. Specifically, 
one proves a functional law of large numbers approximation ("fluid 
limit") for the stochastic processes involved. If the resulting limit- 
ing trajectories satisfy certain properties, then one can conclude that 
the original stochastic system is stable. (For a system described by 
a Markov process, by "stable" we mean positive Harris recurrent.) A 
good exposition can be found in [T7] . This technique is often extended 
by using a Brownian approximation for the underlying stochastic pro- 
cesses of arrivals and service completions. Such assumptions allow a 
description of the queue lengths on a finer ("diffusion") scale, where 
they are non-vanishing; specifically, queues are approximated by a re- 
flected Brownian motion. These techniques present greater technical 
challenges, particularly in deriving steady-state results. We discuss 



diffusion-scaled approximations in Section 3.3, where our exposition 
follows closely the similar results of [TB] . A good reference on diffusion 
approximations of queueing networks is [IB] • 
Next, we examine specific cases. 

3.1. Nj — y oo, (5 fixed. We begin by considering the regime in which 
/3 is fixed, but the arrival rates and the number of chargers tend to 
infinity. More formally, we consider a sequence of systems indexed by 
r, with X\ = rAj and NJ = rNj. (Note that r is a superscript, not 
the r th power.) In this regime, by results from [T3] which we already 
stated above, the average routing rates A[- are guaranteed to stabilize 
the system provided (At) G A, but are not in general guaranteed to 
be socially optimal even as r — y oo. One important result comes from 
considering large systems. If (A») G A, the system is strictly under- 
loaded. Consequently, as r — y oo, the probability of an arrival request 
having to queue and the average queueing time will both converge to 
zero. This is an important consideration in the problem of charging 
electric vehicles. 
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3.2. Nj fixed, — > 0. This is the regime in which the GPD algo- 
rithm converges to the optimal solution of the linear program ([l]). We 
consider a sequence of systems indexed by r, for which j3 r — > and 
all other parameters are unchanged. This corresponds to running a 
sequence of slightly different GPD algorithms on the same external 
pattern of requests. 

All the results for this case, follow straightforwardly from [15j; we 
rewrite them in the notation of our problem. Let Q r At) be the sequence 
of virtual queues. Then, uniformly over any compact time interval of 
t G [0,T], we have 

/3 r Q^(/3 r -) =► qj(-), u.o.c. as (3 r ->■ 0, (2a) 

where the limiting trajectories satisfy 

qj(t) — > q*, as t — > oo. (2b) 

Here, the parameters (q*) are the optimal dual variables corresponding 



to the capacity constraints (lc). 



In addition, we describe the convergence of the routing rates. Let 
X-[j(t) be the exponentially- weighted average of the routing decisions 
up to time t. Then, uniformly on compact sets, 

Xlj(/3 r -) => Xij(-), u.o.c. as /3 r -> 0; (3a) 

where the limiting trajectories satisfy 

Xij(t) — > Xlj, as t — > oo. (3b) 

Here, X*a are the (unique) optimal solution of the linear program ([I]) . 
These results only require the functional law of large numbers given 
in Assumption [TJ but not the functional central limit theorem, as shown 
in |T5]. If the interarrival and service times are assumed to be Poisson, 
and the parameters ^ are rational (so that the virtual queueing system 
is a countable state-space Markov process), these results imply that 
the steady-state routing rates A^ converge to the optimal rates A*,- as 
(5 r -> 0. 

3.3. Nj — > oo, (3 — > 0. In this case, we combine the effects of the pre- 
vious two limits. The operational regime we are interested in is known 
as the Halfin-Whitt regime (THJ. When charging stations have many 
individual chargers, it is possible to operate a heavily loaded system 
while keeping waiting times short and providing a service guarantee on 
the probability of an arriving request having to wait. This effect was 
originally observed in a queueing system with a single pool of many 
servers, but has since been shown to apply in settings with multiple 
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server types, at least under conditions similar to Assumption |4j This 
result leads to the "square- root staffing" principle observed in [19] . 

From the previous sections, we expect to see virtual queues f3Qj ~ q*. 
Our interest will be in the deviations from these values. If (3 is held 
fixed, these deviations may, in general, be large. We will now choose 
— > along an appropriate sequence as the size of the system grows 
larger. We will establish that after rescaling time, the process /3Q(-) can 
be described as a Brownian oscillation around q*. Up to a multiplica- 
tive constant, this result is the best possible, whenever the underlying 
arrival and service completion processes satisfy the functional central 
limit theorem given in Assumption [2] It can be translated into a re- 
flected Brownian motion approximation for the actual queue sizes in 
the system. 

We consider a sequence of systems indexed by r — > oo. The request 
arrival rates satisfy \\ = r\. The sizes of the charger pools satisfy 
iVy = rNj + y/rrij + 0(1), for some collection of values of rij G M. 
(possibly negative); this is the "square-root staffing". (The 0(1) term 
is included to ensure that iVJ is an integer.) We assume that the 
physical limits on the charging rate to not change, so that the charging 
rates \Xij do not change with r. Finally, we choose j3 r = f{r)~ x for 
some function /(r) with r 1 / 2 <C f(r) <C r; for example, we may take 
/(r) = r 3 ' 4 . (Larger values of f(r), corresponding to smaller values of 
(3 r , will lead to more precise but slower convergence.) 

Let Alj (t) be the number of EVs of type % routed to chargers of type 
j during the interval [0, £]. In Theorem [IJ we will establish that this 
quantity can be approximated by A[£ + \/rB(t) for some Brownian 
motion B, provided we consider the system at a time when the GPD 
algorithm has reached its steady-state. The result is similar in nature 
to [IB] , but there are technical differences discussed in the proof. 

We state the next result next, whose proof is given in the Appendix. 

Theorem 1. Assume the arrival process satisfies the FCLT given in 
Assumption^ Let 

q r ,(t)=r^{Q^t)~(PT 1 q*), 

and let 

fiy(*) = r- l/2 (A^t) - A*,i) . 

If(q r j(0)) ->• e 1R J , then d^(-) => H{W), where W is the Brownian 
motion identified in Assumption^ and H is a linear mapping defined in 
Q below. (Thus, H(W) is also a Brownian motion, but with correlated 
components.) 
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Further, suppose the arrival and service completion processes are 
Poisson, and the parameters \iij are rational, so that the virtual queue- 
ing system becomes a countable state-space Markov process. For each 
r consider the associated stationary version of the process. Then, 
(q r j) -^06l J . 

In order to define the load-balancing linear map H, we introduce 
some additional terminology. We call the activities (ij) with Ay > 
basic and let £ be the set of basic activities. By Assumption |4| the 
solution to (fil) is unique. It follows that the bipartite graph with ver- 
tices {car types, charger types} and (undirected) edges corresponding 
to basic activities is acyclic; that is, it is either a tree or a union of 
trees. We can now define the linear map H as follows. For a vector 
v — (vx,...,Vi) G M. 1 , the image w = H(v) (with coordinates indexed 
by basic activities (ij)) satisfies 



£ 



Wij = Vi, Mi 



3 

Ei' w i'i/^i'j __ £i" w i"r 



H'ij H'ij' 



Vi,(ij),(ij'). (4) 



To see that this is a well-defined map, we show how to solve the above 
system of equations. Pick a leaf of one of the connected components of 
the basic activity graph, i.e. a vertex with a single edge coming out of 
it. On this edge, we either have Wij = v^, or can eliminate the variable 
Wij using the second set of equations. Repeating this process, we will 
arrive at a unique solution. 

Next, we discuss the implications of the main result. The conclusion 
of the theorem asserts that the arrival process to each charger type can 
be described as a Brownian oscillation around the optimal rate ^V A*^; 
thus, the arrival process to each charger type satisfies a functional cen- 
tral limit theorem. Because we have assumed (Assumption |2J) that the 
service completion process satisfies it as well, we may approximate the 
number of occupied chargers by a Brownian oscillation about its op- 
timal point, and the queue size as a reflected Brownian motion with 
drift; the drift is given by the collection of parameters nj. In particular, 
the probability of an arriving request being asked to queue will depend 
on the quantity Yl Q.j* n ji the sum being taken over the connected com- 
ponent of the basic activity tree containing the corresponding request 
type. (It will also depend on the linear map H of Ml).) Note that we 
require 0(y/r) overstaffing somewhere in the connected component of 
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each EV type, but do not prescribe where. Consequently, it is advan- 
tageous to arrange the system to have large connected components to 
allow greater freedom in the placement of extra chargers. 



4. Load balancing 

The GPD algorithm will typically produce rates which place a high 
load on one or more of the "least expensive" charger pools. This is 
reasonable in a large system, where square-root overstaffing is suffi- 
cient to deal with the high load, but may be undesirable in a system 
where the charging stations are small. One possibility for mitigating 
this is to reduce the value of Nj in the GPD algorithm, so that on 
average only a fraction of the chargers may be used. However, this 
reduces the stability region of the algorithm. A better way to mitigate 
this difficulty is to encourage the charging stations to spread the load 
more evenly. Group the charging stations into clusters. (Clusters may 
overlap.) Change the objective function of Q to 

minimize ^ \jCi(j) + ^ Wipi (5a) 

s.t. Xi = yjAy, Vz (5b) 



ty > E — * ^ ( 5c ) 

over Xij > 0, Vi, j. (5e) 

where I runs over the clusters, Wi are weights, and pi is the maximal 
load of any charging station in the cluster. Modifying the objective in 
this manner means that charging stations within each cluster will "try" 
to have equal loads, because only the maximal load within a cluster 
is penalized. By adjusting the weight vector Wi, we can change the 
relative importance of spreading the load across different stations, and 
finding the lowest-cost routing pattern. 

The corresponding modification of the GPD algorithm is called the 
Load Balancing (LB) algorithm. 



LB Algorithm: 
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1. Each charger type maintains a virtual queue variable Qj(t), 
initialized e.g. to Qj(t) = 0. Each charger type also maintains 
a virtual queue Lj(t) for each cluster L to which j belongs. 
These are initialized so that ^2j & iLj(t) = Wi for each /. (We 
assume that this is possible.) 

2. When a vehicle of type i requests service at time t, we locate 



j* G argmincj(j) + 



PiQAt) + L&)) 



J fJ'ij 

The vehicle announces its decision to j*, and the corresponding 
virtual queues are incremented: 

Q* H- Q* + —, L* ^ L* + — . 

fMj* r'ij* 

3. All virtual queues Qj are decreased at a rate Nj per time unit 
whenever they are positive. The cluster virtual queues Lj for 
j G I are decreased when Yljei 0Lj > Wj; if this is the case, we 
decrease 

Lj i-> Lj - Nj, Vj G /. 

Note that the algorithm is still distributed, since additional communi- 
cation needs to happen only between the charging stations within the 
cluster. Namely, stations must communicate their updated values of 
Lj whenever these increase, and some entity must decide to decrease 
the Lj within the cluster once their sum exceeds j3~ 1 Wi. 

Similarly to the analysis of GPD, as /3 — > the routing pattern 
produced by the LB algorithm will converge to the optimal solution to 
d5l). The virtual queues (3Lj will converge to the optimal dual variables 



corresponding to the constraints defining the loads, (5d). Further, 
under Assumptions [2] and |4j a variant of Theorem [T] can be shown 
to hold, with deviations of the routing process from its nominal value 
given by some Brownian motion; however, the Brownian motion will 
now be more cumbersome to define. 

In Section [6j numerical work shows that replacing the GPD algo- 
rithm by the LB algorithm can lead to substantial reduction in the 
queueing delays encountered by the vehicles. 

A feature of the GPD and LB algorithms is that they can be slow 
to reach convergence. Specifically, the typical time scale of GPD and 
LB is /3 _1 ; so as /3 — > 0, the algorithms get more precise but slower at 
reacting to shocks. A possibility to mitigate this is to run GPD (or LB) 
in the background to identify the set of edges in the basic activity tree, 
but then run a different algorithm once that tree has been identified. 
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We will return to this idea in Section [5j where we present one such 
possible algorithm. 



5. The case of No costs 

We now consider the special case of Cj(j) 6 {0, oo} for all i and j, 
corresponding to vehicles that are indifferent to the choice of charging 
station provided they can reach it (and it has compatible technology). 
This corresponds to the question of optimizing the throughput, i.e. the 
number of vehicles that are served. In this setting, the parameter /3 in 
the GPD algorithm is irrelevant, because the decisions will be the same 
for all values of j3. In particular, we conclude that the GPD algorithm 
always generates routing rates in the feasible region, provided A is 
nonempty. 

Note that if the feasible region is nonempty, the optimal dual vari- 



ables corresponding to the capacity constraints (lc) are q* = 0. While 
of course we will not have Qj (t) = at all times for the virtual queues 
under GPD, they will be close to 0, and will periodically hit unless 
the system is overloaded. 

If the interarrival and service times are exponential, there will be 
well-defined (unique) steady-state routing rates A*,-, which we know 
are feasible for ([I]). However, if Ci(j) E {0, oo}, the optimal solution to 
(fil) will certainly not be unique. Determining to which optimal solution 
the routing rates of GPD converge is difficult; this is another reason to 
use some form of load-balancing. 

Note that for the GPD algorithm, in the case of zero costs, conver- 
gence time is not an issue: routing rates are "always" feasible. However, 
for the LB one, the quantity (3 makes an appearance, since we decre- 
ment the cluster queues Lj when Yljei @Lj = Wi. Picking values of 
/3 that are too small will result in imperfect load balancing (although 
typically will reduce the maximal loads somewhat). A better, and 
faster, possibility is to use the Freest Charger Shortest Queue (FCSQ) 
algorithm, modelled on [20] . 



FCSQ Algorithm: 

1. Identify the basic activity tree (e.g. by running GPD in the 
background) . 

2. When a vehicle of type % requests service, if some charger with 
Ci(j) < oo and //^ > is available, route the vehicle to the 
charging station with the largest fraction of free chargers. 
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3. If no charger is available in a charging station, the charging 
station sends back an estimate of the time when the vehicle will 
enter queueing. The vehicle joins the charging station where 
this time is earliest. 

The last step of the algorithm is deliberately vague, because many 
versions should give essentially similar results. In the simulation pre- 
sented in Section |6j we add up the (future) charging times of all of the 
EVs in the queue and divide by the number of chargers at the station. 
We could add to that the residual charging times of the vehicles that 
already charging; or we could replace the exact times by estimates, re- 
porting iV~ 1 (^//^ m), where n, is the number of queued vehicles of 
type i. 

The FCSQ algorithm should guarantee that, when there are free 
chargers, they are divided fairly among the various charger pools, while 
when there are queues, they have similar waiting times at all the charg- 
ers. This means that when there is a sharp spike in arrivals at one of 
the charger pools, these arrivals are quickly spread out as far as pos- 
sible, helping the spike dissipate faster. Strong analytic guarantees 
are available for a similar algorithm (Longest Queue Freest Server, see 
|20j ) in the case when the non-zero charging rates pij depend only on 
the charging station technology (i.e. pij = /ij), and partial results are 
available for general parameters. Similar algorithms are common in the 
queueing literature, see for example |21j . or [22] and references therein. 



6. Numerical Illustration 

We now consider the toy network shown in Fig. [I] to illustrate the 
performance of the various algorithms presented in this work. There 
are two EV classes, three charger classes, and each EV class can use 
two of the charger classes. The service rates \Xij are 

[1 3 
fX ~ [o 1 2 

We use Nj = 20 for all j. We simulate 10,000 EV arrivals. The first 
5,000 arrivals are generated using arrival rates A = (2.5,2.2) x 20; 
for the second 5,000 arrivals, we reverse the arrival rates to obtain 
A = (2.2, 2.5) x 20. We do this to illustrate the effect of a change 
in the arrival pattern that does not change the basic activity tree on 
the algorithms under consideration. Note that the system is heavily 
loaded, but not overloaded: the load-balancing solution to the linear 
program ([5]), achieved for sufficiently high values of W\, has p = 0.91 
for the first arrival pattern, and p = 0.97 for the second arrival pattern. 
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did 



20 20 20 

Figure 1. The diagram for the toy example indicating 
Nj and /%. 

In a conventional single-server queue, such heavy loads would result in 
long queues and very high probability of having to queue; in the many- 
charger setting, however, approximately 70% of the arriving vehicles 
are taken into service without having to queue at all. 

Fig.[2]presents the delays observed when running the GPD algorithm. 
The three curves present the delays encountered by the vehicles routed 
to each of the three charging stations. (One of the curves stays at 
throughout the simulation.) The vertical line indicates the time when 
the arrival pattern changes. We see that before the change in the 
arrival pattern, a large queue had been forming at station 2 (dashed 
line); after the change, it slowly goes away, but a sizeable queue forms 
at station 3 (dotted line) instead. 

In Fig. [3j we simulate the LB algorithm, putting all of the chargers 
into a single cluster. We use W\/ (3 = 1000. The scale is the same as for 
the GPD routing; so we see that the delays have gotten shorter. (The 
maximal delay has dropped from 1.64 to 1.3.) We also see the marked 
change in the pattern after the change in the arrival rate pattern: the 
delay at station 2 (dashed line) drops nearly to zero, however, it takes 
a while for that to occur. The behavior of the LB algorithm would be 
improved if we considered a larger system (Nj > 20). 

Fig.[4]shows the delays arisint from FCSQ routing on the same arrival 
data. Note that the largest delay is now only 0.94, more than a third 
smaller than under GPD. Also, when there is queueing delay, it is 
nearly the same at all stations, meaning that there is little incentive 
for any one vehicle to disobey the algorithm and go elsewhere. This 
is in shaprp contrast to GPD, which tends to produce highly unequal 
queueing delays. 

7. Discussion and future work 

We have constructed a set of algorithms - GPD, LB, and FCSQ - 
which collectively route the EVs to chargers in the network in such 
a way as to avoid long queueing and adjust to changing demand pat- 
terns whenever this is possible. The LB algorithm s designed to balance 
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Figure 2. Delays produced by GPD routing 



the loads across different charging stations; the FCSQ algorithm cor- 
responds to a faster version of LB, which, however, can only be run 
after LB or GPD has been used to identify the basic activities along 
which routing should happen. All of the algorithms are distributed, 
and therefore can scale to a large network. 

The algorithms have implications for building the charging network. 



As we observed in Section |3.3[ square-root overstaffing is needed for 
the operation of a system with many chargers in each pool and a small 
value of the parameter /3 with short queueing delays. Moreover, the 
overstaffing is aggregated over all of the charging stations in a single 
connected component of the basic activity graph. Therefore, for greater 
freedom in the placement of extra charging facilities it is desirable to 
keep these connected components large. That is, in contrast to the 
facility locations problem, we not only allow, but encourage EVs of a 
single type to visit multiple charger types. 

Two important extensions of this work are of interest. First, it is 
desirable to arrange charging stations so that the basic activity graph 
for the system will have large connected components. However, this 
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Figure 3. Delays produced by LB routing 



is a highly nonlinear and nonconvex requirement. We would like to 
create an algorithm for generating arrangements of charging stations 
satisfying this condition. Second, we would like to incorporate the ef- 
fect of batteries onto this system. Because of the high load on the 
electrical network produced by charging an electric vehicle, it is desir- 
able to add battery capacity to the charging stations to smooth the 
peak demand. This introduces an additional algorithmic challenge by 
introducing demand mobility in time in addition to in space. 

Finally, we remark that due to technical challenges, some of the 
results concerning the Halfm-Whitt regime remain conjectural. Specif- 
ically, there is a gap between the finite-time-horizon results on the 
Brownian motion approximations, and steady-state quantities such as 
the probability of an arriving request being asked to queue. However, 
this gap is not very important in practice because of the diurnal de- 
mand pattern. 

[Proof of Theorem [I] The overall technique of proof of Theorem [I] 
follows the proof of [16, Theorem 6.3]. We outline both the argument 
and the necessary changes below. 
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Figure 4. Delays produced by FCSQ routing 



Sketch of proof of Theorem [1} To show the first assertion of the theo- 
rem, we begin by demonstrating that if the virtual queues are initially 
close to their equilibrium values, then the deviations must remain small 
throughout the entire time interval [0, i\. The version of this argument 
in [THl Theorem 6.1] does not hold in our case; however, the following 
result is true. 



i„* 



Lemma 2 (Version of Theorem 6.1 of [16]). If QJ(0) - (/3TV 
o(^), then uniformly on any compact set t G [0, T] we have 

(3 r Q V j(t) — q* = o(l), Vj, uniformly on any compact set t G [0,T]. 

Consequently, 



CiO) + — - — - c ih ) + 



r'ij J V r'ij' 

uniformly on any compact set t G [0, T]. 



o(y/r), V(ij),(ij')eE, 



We postpone the proof of the lemma to proceed with the rest of the 
argument. The first claim of the lemma implies that, for r sufficiently 
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large, vehicles will not be routed along non-basic activities, because dif- 
ferences in virtual queue sizes dominate the fluctuations in the size of 
each virtual queue. The second claim then determines the proportion 
of vehicles routed along each of the basic activities. Because the differ- 
ences between the quantities controlling the routing are much smaller 
than the size of the Brownian fluctuations in the arrivals, the Brownian 
fluctuations are unaffected by them. 

The stationary result is proved by considering the fluid limits of the 
form \/r(Q r (y/ru) — q])- The analysis there shows that in steady state, 
this deviation is o(y/r); indeed, it should be o((/3' r ) _1 ). We refer the 
reader to [TBI Theorem 6.3, 6.4] for details. □ 

We now outline the proof of Lemma [2} 

Sketch of proof of Lemma [£| Both claims of the lemma are proved us- 
ing the technique of local fluid limits. For details, the reader may 
consult e.g. [221 Section 8]. 

For the first assertion, define for all j the local-fluid-scaled queueing 
process 

qf' m) (u)=PQ r j (f n + Pu)-q* j , UE[0,T], m = l,...,r/7. 

It is standard to show that, as r — y oo, each sequence q f (•) con- 
verges to the family of Lipschitz processes q™(-) satisfying certain dif- 
ferential equations whenever they are defined. Our goal is to show that 
max, qj (u) decreases whenever it is positive, and does so at some 
rate bounded away from 0. Note first that for finite values of g( r,m ), 
the values of Q r are still almost proportional to their nominal val- 
ues q*, and therefore all requests must be routed along basic activities 
only. Moreover, consider the set of charger types that have the largest 
value of q m (u). The routing rule of GPD algorithm ensures that these 
charger types will have only those vehicles routed to them that cannot 
be routed elsewhere. By Assumption [4], if the maximal value of q m 
is positive, then the virtual queues for which q m is maximal will have 
arrivals that are smaller than nominal, and of course will be decreased 
at the nominal rate. Consequently, when the maximal value of q m is 
positive, it must decrease towards its equilibrium value of 0. The time 
T is picked to be sufficiently large so that, e.g., if q™{$) < 1 for all 
j, then q^iT) = for all j. We obtain that on a single time interval, 
q( r ' m '(u) is very likely to stay close to 0. Finally, we show that this 
holds on all r(3 r subintervals simultaneously. Suppose not, and pick 
the subsequence of intervals m(r) along which g( r ' m ) first crosses level 
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e < 1. Pick a subsequence along which q(' r ' m \0) converges, and ob- 
serve that the corresponding local fluid limit contradicts our previous 
assertions. 

For the second assertion, we define the slightly different local-fluid- 
scaled processes 

q( r ' m \ u )= r - 1 / 2 (/3 r Q r J (t m + r- 1 / 2 u)-q*), u E [0, T], m = 1, . . . , y/f. 

Very similar arguments tell us that for finite values of q( r ' m \ the quan- 
tities Ci(j) + Hij 1 (3jQ r j are approximately proportional to their nominal 
values, and therefore routing decisions are very close to nominal. This 
easily implies that the quantities q(j) + jj^ 1 /3jQ r j stay close to their 
nominal values. □ 
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